SETUP

Installing the following packages if not already present

LOADING LIBRARIES

Loading the required libraries

Read in Phenotype file with SRR numbers

pheno_file =  read.table("results/GSE95632_Phenotype_withoutQC.txt", sep="\t")
colnames(pheno_file) <- c("SRR","SampleID","ID","Species","Tissue","Treatment","genome","Antibody","Ab+Treatment")
DT::datatable(pheno_file, rownames=FALSE, options = list(columnDefs = list(list(className = 'dt-center', targets = "_all"))))

Subset samples and add other columns to make .csv file

#Remove input controls
pheno.main <- pheno_file[-c(2,5,17,18),] %>%select("SampleID")

#Peak file paths
Ab = as.factor(gsub(" .*","",pheno.main$SampleID))
Code = as.factor(gsub(".*_","",pheno.main$SampleID))
s_names <- as.factor(paste(Ab,Code,sep = "_"))

#Wrangle and edit columns to get necessary format
#SampleID
pheno.main$SampleID <- as.factor(gsub(".*_","",pheno.main$SampleID))
#Tissue
pheno.main$Tissue <- as.factor("ASM")
#Factor
pheno.main$Factor<- as.factor(pheno_file$Antibody[-c(2,5,17,18)])     
#Treatment
pheno.main$Condition <- as.factor(c("control","control","dex","dex","dex_GFP","dex_GFP","dex_KLF","dex_KLF","control","control","dex","dex","dex_GFP","dex_GFP","dex_KLF","dex_KLF"))
pheno.main$Treatment <- as.factor(pheno.main$Condition)
#Replicates
pheno.main$Replicate <- as.factor(c(1,2))
#ControlID
pheno.main$ControlID <- as.factor(c("1D_input","1D_input","2A_input","2A_input","2A_input","2A_input","2A_input","2A_input"))

#File paths
srr_nos <- pheno_file[-c(2,5,17,18),]%>%select("SRR")
pheno.main$bamReads <- as.factor(paste("/Users/diwadkar/Desktop/ChIP-Seq/GSE95632/data/bam/",srr_nos$SRR,".bam",sep=""))
pheno.main$bamControl <- as.factor(paste("/Users/diwadkar/Desktop/ChIP-Seq/GSE95632/data/bam/",c("SRR5309352","SRR5309352","SRR5309355","SRR5309355","SRR5309355","SRR5309355","SRR5309355","SRR5309355"),".bam",sep=""))
pheno.main$Peaks <- as.factor(paste("/Users/diwadkar/Desktop/ChIP-Seq/GSE95632/data/macs2/",s_names,"_peaks.narrowPeak",sep=""))
pheno.main$PeakCaller <- as.factor("bed")
rownames(pheno.main) <- seq(length=nrow(pheno.main)) 

#View
#GR
pheno.gr <- pheno.main %>% filter(Factor=="GR")
DT::datatable(pheno.gr, rownames=FALSE, options = list(columnDefs = list(list(className = 'dt-center', targets = "_all"))))
write.csv(pheno.gr,"~/Desktop/ChIP-Seq/GSE95632/GSE95632info_gr.csv",row.names=FALSE)
#RNAP2
pheno.rp <- pheno.main %>% filter(Factor=="RNAP2")
DT::datatable(pheno.rp, rownames=FALSE, options = list(columnDefs = list(list(className = 'dt-center', targets = "_all"))))
write.csv(pheno.rp,"~/Desktop/ChIP-Seq/GSE95632/GSE95632info_rp.csv",row.names=FALSE)

Obtaining differentially bound sites

Reading in the peaksets

Set up an experiment to analyze is with a sample sheet (.csv). The peaksets are read in using the dba function.

The metadata shows how many peaks (108812 sites) after merging overlapping ones (140789 total), and the dimensions of dba.plotPCA (default binding matrix of 16 samples by the N sites that overlap in at least two of the samples).

A correlation heatmap gives an initial clustering of the samples using the cross-correlations of each row of the binding matrix.

#GR
dat.gr <- dba(sampleSheet="/Users/diwadkar/Desktop/ChIP-Seq/GSE95632/GSE95632info_gr.csv") 
## 1D ASM GR control control 1 bed
## 1E ASM GR control control 2 bed
## 2A ASM GR dex dex 1 bed
## 2E ASM GR dex dex 2 bed
## 3D ASM GR dex_GFP dex_GFP 1 bed
## 3E ASM GR dex_GFP dex_GFP 2 bed
## 4B ASM GR dex_KLF dex_KLF 1 bed
## 4C ASM GR dex_KLF dex_KLF 2 bed
plot(dat.gr)

#RNAP2
dat.rp <- dba(sampleSheet="/Users/diwadkar/Desktop/ChIP-Seq/GSE95632/GSE95632info_rp.csv") 
## 5E ASM RNAP2 control control 1 bed
## 5b ASM RNAP2 control control 2 bed
## 6d ASM RNAP2 dex dex 1 bed
## 6e ASM RNAP2 dex dex 2 bed
## 7c ASM RNAP2 dex_GFP dex_GFP 1 bed
## 7e ASM RNAP2 dex_GFP dex_GFP 2 bed
## 9b ASM RNAP2 dex_KLF dex_KLF 1 bed
## 9e ASM RNAP2 dex_KLF dex_KLF 2 bed
plot(dat.rp)

Counting reads

The next step is to calculate a binding matrix with scores based on read counts for every sample (affinity scores), rather than confidence scores for only those peaks called in a specific sample (occupancy scores).

Use dba.count function. The current study is based on a transcription factor GR that binds to the DNA, resulting in “punctate”, narrow peaks, it is advisable to use the “summits” option to re-center each peak around the point of greatest enrichment. This keeps the peaks at a consistent width (in this case, with summits=250, the peaks will be 500bp, extending 250bp up- and down- stream of the summit).

If you do not have the raw reads available to you, this object is available loading using e.g. data(tamoxifen_counts).

#GR
dat.gr <- dba.count(dat.gr,summits=250) # This step takes a while
## Re-centering peaks...
#RNAP2
dat.rp <- dba.count(dat.rp,summits=250)
## Re-centering peaks...

This shows that all the samples are using the same, 82032 length consensus peakset. Also, a new column has been added, called FRiP (Fraction of Reads in Peaks). This is the proportion of reads for that sample that overlap a peak in the consensus peakset, and can be used to indicate which samples show more enrichment overall.

dat.gr
## 8 Samples, 98680 sites in matrix:
##   ID Tissue Factor Condition Treatment Replicate Caller Intervals FRiP
## 1 1D    ASM     GR   control   control         1 counts     98680 0.30
## 2 1E    ASM     GR   control   control         2 counts     98680 0.30
## 3 2A    ASM     GR       dex       dex         1 counts     98680 0.31
## 4 2E    ASM     GR       dex       dex         2 counts     98680 0.35
## 5 3D    ASM     GR   dex_GFP   dex_GFP         1 counts     98680 0.29
## 6 3E    ASM     GR   dex_GFP   dex_GFP         2 counts     98680 0.28
## 7 4B    ASM     GR   dex_KLF   dex_KLF         1 counts     98680 0.28
## 8 4C    ASM     GR   dex_KLF   dex_KLF         2 counts     98680 0.25
dat.rp
## 8 Samples, 29852 sites in matrix:
##   ID Tissue Factor Condition Treatment Replicate Caller Intervals FRiP
## 1 5E    ASM  RNAP2   control   control         1 counts     29852 0.09
## 2 5b    ASM  RNAP2   control   control         2 counts     29852 0.10
## 3 6d    ASM  RNAP2       dex       dex         1 counts     29852 0.11
## 4 6e    ASM  RNAP2       dex       dex         2 counts     29852 0.10
## 5 7c    ASM  RNAP2   dex_GFP   dex_GFP         1 counts     29852 0.09
## 6 7e    ASM  RNAP2   dex_GFP   dex_GFP         2 counts     29852 0.06
## 7 9b    ASM  RNAP2   dex_KLF   dex_KLF         1 counts     29852 0.08
## 8 9e    ASM  RNAP2   dex_KLF   dex_KLF         2 counts     29852 0.09

Plot a new correlation heatmap based on the affinity scores.

plot(dat.gr)

plot(dat.rp)

Establishing a contrast

Before running the differential analysis, we need to tell DiffBind which cell lines fall in which groups.

This uses the Condition metadata (dex vs. control) to set up a contrast with 2 samples in the dex group and 2 samples in control group.

The comparison indicates the group 1 (dex) vs. group 2 (control)

#GR
dat.gr <- dba.contrast(dat.gr, group1=dat.gr$masks$dex, group2=dat.gr$masks$control, name1="dex", name2="control")
dat.gr <- dba.contrast(dat.gr, group1=dat.gr$masks$dex_GFP, group2=dat.gr$masks$dex_KLF, name1="dex_GFP", name2="dex_KLF")
#dat <- dba.contrast(dat, categories=DBA_CONDITION, minMembers=2)

#RNAP2
dat.rp <- dba.contrast(dat.rp, group1=dat.rp$masks$dex, group2=dat.rp$masks$control, name1="dex", name2="control")
dat.rp <- dba.contrast(dat.rp, group1=dat.rp$masks$dex_GFP, group2=dat.rp$masks$dex_KLF, name1="dex_GFP", name2="dex_KLF")
#dat <- dba.contrast(dat, categories=DBA_CONDITION, minMembers=2)

Performing the differential analysis

The default analysis will run an DESeq2 analysis using the default binding matrix. The resultant DBA object shows that 37582 of the 82032 sites are identified as being significantly differentially bound (DB) using a FDR <= 0.05. A correlation heatmap can be plotted, based on the result of the analysis.

#GR - dex vs control
dat.gr <- dba.analyze(dat.gr)
dat.gr
## 8 Samples, 98680 sites in matrix:
##   ID Tissue Factor Condition Treatment Replicate Caller Intervals FRiP
## 1 1D    ASM     GR   control   control         1 counts     98680 0.30
## 2 1E    ASM     GR   control   control         2 counts     98680 0.30
## 3 2A    ASM     GR       dex       dex         1 counts     98680 0.31
## 4 2E    ASM     GR       dex       dex         2 counts     98680 0.35
## 5 3D    ASM     GR   dex_GFP   dex_GFP         1 counts     98680 0.29
## 6 3E    ASM     GR   dex_GFP   dex_GFP         2 counts     98680 0.28
## 7 4B    ASM     GR   dex_KLF   dex_KLF         1 counts     98680 0.28
## 8 4C    ASM     GR   dex_KLF   dex_KLF         2 counts     98680 0.25
## 
## 2 Contrasts:
##    Group1 Members1  Group2 Members2 DB.DESeq2
## 1     dex        2 control        2     39099
## 2 dex_GFP        2 dex_KLF        2     28861
plot(dat.gr, contrast=1)

plot(dat.gr, contrast=2)

#GR - GFP dex vs KLF dex
#dat.gr2 <- dba.analyze(dat.gr2)
#dat.gr2
#plot(dat.gr2, contrast=1)

#RNAP2 - dex vs control
dat.rp <- dba.analyze(dat.rp)
dat.rp
## 8 Samples, 29852 sites in matrix:
##   ID Tissue Factor Condition Treatment Replicate Caller Intervals FRiP
## 1 5E    ASM  RNAP2   control   control         1 counts     29852 0.09
## 2 5b    ASM  RNAP2   control   control         2 counts     29852 0.10
## 3 6d    ASM  RNAP2       dex       dex         1 counts     29852 0.11
## 4 6e    ASM  RNAP2       dex       dex         2 counts     29852 0.10
## 5 7c    ASM  RNAP2   dex_GFP   dex_GFP         1 counts     29852 0.09
## 6 7e    ASM  RNAP2   dex_GFP   dex_GFP         2 counts     29852 0.06
## 7 9b    ASM  RNAP2   dex_KLF   dex_KLF         1 counts     29852 0.08
## 8 9e    ASM  RNAP2   dex_KLF   dex_KLF         2 counts     29852 0.09
## 
## 2 Contrasts:
##    Group1 Members1  Group2 Members2 DB.DESeq2
## 1     dex        2 control        2      1616
## 2 dex_GFP        2 dex_KLF        2       235
plot(dat.rp, contrast=1)

plot(dat.rp, contrast=2)

#RNAP2 - GFP dex vs KLF dex
#dat.rp2 <- dba.analyze(dat.rp2)
#dat.rp2
#plot(dat.rp2, contrast=1)

Retrieve the differentially bound sites

The output is a GRanges object, appropriate for downstream processing.

The value columns show:

  • the mean read concentration over all the samples: the default calculation uses log2 normalized ChIP read counts with control read counts subtracted
  • the mean concentration of the first (dex) group
  • the mean concentration of second (control) group
  • the Fold column shows the difference in mean concentrations between the two groups (Conc_dex - Conc_control): a positive value indicating increased binding affinity in the dex group and a negative value indicating increased binding affinity in the control group.
  • a raw p-value: confidence measures for identifying these sites as differentially bound
  • FDR: a multiple testing corrected FDR in the final column.
#GR1
dat.gr.DB <-  dba.report(dat.gr)
dat.gr.DB
## GRanges object with 39099 ranges and 6 metadata columns:
##         seqnames              ranges strand |      Conc  Conc_dex
##            <Rle>           <IRanges>  <Rle> | <numeric> <numeric>
##   17705    chr11 114179230-114179730      * |      9.44     10.41
##   11837    chr10   74056322-74056822      * |      9.12     10.06
##   42623     chr2   20438541-20439041      * |      8.75      9.67
##   17699    chr11 114154432-114154932      * |      8.78      9.72
##   36799    chr17   59825482-59825982      * |      8.69      9.61
##     ...      ...                 ...    ... .       ...       ...
##    8163     chr1 221989640-221990140      * |      5.86       5.4
##   73147     chr5 137188108-137188608      * |      2.83      3.61
##   88784     chr8   67245151-67245651      * |      2.83      3.61
##   11693    chr10   71230593-71231093      * |       6.1      6.41
##   25151    chr13   74229142-74229642      * |      3.35      4.04
##         Conc_control      Fold   p-value       FDR
##            <numeric> <numeric> <numeric> <numeric>
##   17705         4.71       5.7 1.16e-119 1.15e-114
##   11837         5.36       4.7 2.93e-111 1.44e-106
##   42623         5.38      4.29  6.53e-93  2.15e-88
##   17699         5.26      4.46  4.36e-92  1.08e-87
##   36799         5.51       4.1  7.36e-83  1.45e-78
##     ...          ...       ...       ...       ...
##    8163          6.2      -0.8    0.0198      0.05
##   73147         0.98      2.63    0.0198      0.05
##   88784         0.98      2.63    0.0198      0.05
##   11693         5.69      0.72    0.0198      0.05
##   25151         1.96      2.08    0.0198      0.05
##   -------
##   seqinfo: 33 sequences from an unspecified genome; no seqlengths
#GR2
#dat.gr2.DB <-  dba.report(dat.gr2)
#dat.gr2.DB

#RNAP2
dat.rp.DB <-  dba.report(dat.rp)
dat.rp.DB
## GRanges object with 1616 ranges and 6 metadata columns:
##         seqnames              ranges strand |      Conc  Conc_dex
##            <Rle>           <IRanges>  <Rle> | <numeric> <numeric>
##   17681    chr20   47318285-47318785      * |      7.84      8.63
##   10608    chr16   56608395-56608895      * |       9.2      9.79
##   18463    chr22   30246411-30246911      * |      7.36      5.79
##   23916     chr6   35727765-35728265      * |      6.22       7.1
##   28928     chr9 129320831-129321331      * |      7.59      8.26
##     ...      ...                 ...    ... .       ...       ...
##   23416     chr5 179829909-179830409      * |      4.26      3.07
##   10695    chr16   66996370-66996870      * |      3.43      4.23
##    3982    chr10 110535930-110536430      * |      3.41      1.32
##   17476    chr20   35443015-35443515      * |      3.59       4.4
##   20292     chr3 155079695-155080195      * |      4.16      4.83
##         Conc_control      Fold   p-value       FDR
##            <numeric> <numeric> <numeric> <numeric>
##   17681         5.98      2.64  1.95e-36  5.47e-32
##   10608         8.17      1.63  2.52e-25  3.55e-21
##   18463          8.1      -2.3  4.13e-24  3.87e-20
##   23916         3.65      3.45  1.94e-19  1.36e-15
##   28928         6.29      1.97  2.59e-18  1.46e-14
##     ...          ...       ...       ...       ...
##   23416          4.9     -1.83   0.00286    0.0499
##   10695         1.43       2.8   0.00286    0.0499
##    3982         4.23     -2.91   0.00287    0.0499
##   17476         1.53      2.87   0.00287    0.0499
##   20292         2.84      1.99   0.00287    0.0499
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths
#RNAP2
#dat.rp2.DB <-  dba.report(dat.rp2)
#dat.rp2.DB

Visualization

PCA

A PCA plot shows a better understanding of clustering and association between samples than heatmaps.

#PCA plot for normalized read counts for all binding sites
dba.plotPCA(dat.gr,DBA_TREATMENT,label=DBA_CONDITION)

#dba.plotPCA(dat.gr2,DBA_TREATMENT,label=DBA_CONDITION)
#dba.plotPCA(dat,DBA_FACTOR,label=DBA_CONDITION)
#dba.plotPCA(dat.rp1,DBA_TREATMENT,label=DBA_CONDITION)
dba.plotPCA(dat.rp,DBA_TREATMENT,label=DBA_CONDITION)

#PCA plot using only the differentially bound sites, using an FDR threshold of 0.05
dba.plotPCA(dat.gr,contrast = 1,label=DBA_REPLICATE)

dba.plotPCA(dat.gr,contrast = 2,label=DBA_REPLICATE)

#dba.plotPCA(dat.gr1,contrast = 1,label=DBA_TREATMENT)
#dba.plotPCA(dat.gr2,contrast = 1,label=DBA_REPLICATE)
#dba.plotPCA(dat.gr2,contrast = 1,label=DBA_TREATMENT)
#dba.plotPCA(dat,contrast = 1,label=DBA_FACTOR)

#RNAP2
dba.plotPCA(dat.rp,contrast = 1,label=DBA_REPLICATE)

dba.plotPCA(dat.rp,contrast = 2,label=DBA_REPLICATE)

#dba.plotPCA(dat.rp1,contrast = 1,label=DBA_TREATMENT)
#dba.plotPCA(dat.rp2,contrast = 1,label=DBA_REPLICATE)
#dba.plotPCA(dat.rp2,contrast = 1,label=DBA_TREATMENT)

MA PLOTS

MA plots are a useful way to visualize the effect of normalization on data, as well as seeing which of the datapoints are being identified as differentially bound. An MA plot can be obtained for the control vs dex treatment cell:

dba.plotMA(dat.gr)

dba.plotMA(dat.rp)

#Same data shown with the concentrations of each sample groups plotted against each other - 
dba.plotMA(dat.gr, bXY=TRUE)

#dba.plotMA(dat.gr2, bXY=TRUE)
dba.plotMA(dat.rp, bXY=TRUE)

#dba.plotMA(dat.rp2, bXY=TRUE)

VOLCANO PLOTS

Volcano plots also highlight significantly differentially bound sites and show their fold changes.

dba.plotVolcano(dat.gr)

#dba.plotVolcano(dat.gr2)
dba.plotVolcano(dat.rp)

#dba.plotVolcano(dat.rp2)

BOXPLOTS

Box plots of read distributions for significantly differentially bound (DB) sites - The plot shows in the first two boxes that amongst differentially bound sites overall, the dex samples have a somewhat higher mean read concentration. The next two boxes show the distribution of reads in differentially bound sites that exhibit increased affinity in the dex samples, while the final two boxes show the distribution of reads in differentially bound sites that exhibit increased affinity in the Control samples.

Pvals is a matrix of p-values (computed using a two-sided Wilcoxon ‘MannWhitney’ test, paired where appropriate) indicating which of these distributions are significantly different from another distribution.

#GR1
sum(dat.gr.DB$Fold<0)
## [1] 15983
sum(dat.gr.DB$Fold>0)
## [1] 23116
pvals1 <- dba.plotBox(dat.gr)

pvals1
##                dex.DB control.DB   dex.DB+ control.DB+   dex.DB-
## dex.DB       1.00e+00   0.00e+00 2.82e-212           0 9.73e-131
## control.DB   0.00e+00   1.00e+00 4.24e-120           0  0.00e+00
## dex.DB+     2.82e-212  4.24e-120  1.00e+00           0  0.00e+00
## control.DB+  0.00e+00   0.00e+00  0.00e+00           1  0.00e+00
## dex.DB-     9.73e-131   0.00e+00  0.00e+00           0  1.00e+00
## control.DB-  0.00e+00   0.00e+00  0.00e+00           0  0.00e+00
##             control.DB-
## dex.DB                0
## control.DB            0
## dex.DB+               0
## control.DB+           0
## dex.DB-               0
## control.DB-           1
#GR2
#sum(dat.gr2.DB$Fold<0)
#sum(dat.gr2.DB$Fold>0)

#pvals2 <- dba.plotBox(dat.gr2)
#pvals2

#RNAP2-1
sum(dat.rp.DB$Fold<0)
## [1] 576
sum(dat.rp.DB$Fold>0)
## [1] 1040
pvals3 <- dba.plotBox(dat.rp)

pvals3
##               dex.DB control.DB   dex.DB+ control.DB+   dex.DB-
## dex.DB      1.00e+00   2.58e-33  1.21e-74    7.23e-22  4.39e-35
## control.DB  2.58e-33   1.00e+00  4.67e-17    1.86e-52  2.65e-83
## dex.DB+     1.21e-74   4.67e-17  1.00e+00    4.96e-96 3.85e-155
## control.DB+ 7.23e-22   1.86e-52  4.96e-96    1.00e+00  4.97e-01
## dex.DB-     4.39e-35   2.65e-83 3.85e-155    4.97e-01  1.00e+00
## control.DB- 2.12e-90   5.44e-25  1.41e-02   1.18e-108 9.99e-172
##             control.DB-
## dex.DB         2.12e-90
## control.DB     5.44e-25
## dex.DB+        1.41e-02
## control.DB+   1.18e-108
## dex.DB-       9.99e-172
## control.DB-    1.00e+00
#RNAP2-2
#sum(dat.rp2.DB$Fold<0)
#sum(dat.rp2.DB$Fold>0)

#pvals4 <- dba.plotBox(dat.rp2)
#pvals4

HEATMAPS

#GR
#Simple heatmap
corvals1 <- dba.plotHeatmap(dat.gr)

#Binding affinity heatmap
corvals1 <- dba.plotHeatmap(dat.gr, contrast=1, correlations=FALSE)

corvals1 <- dba.plotHeatmap(dat.gr, contrast=2, correlations=FALSE)

#Simple heatmap
#corvals2 <- dba.plotHeatmap(dat.gr2)

#Binding affinity heatmap
#corvals2 <- dba.plotHeatmap(dat.gr2, contrast=1, correlations=FALSE)

#RNAP2
#Simple heatmap
corvals3 <- dba.plotHeatmap(dat.rp)

#Binding affinity heatmap
corvals3 <- dba.plotHeatmap(dat.rp, contrast=1, correlations=FALSE)

corvals3 <- dba.plotHeatmap(dat.rp, contrast=2, correlations=FALSE)

#Simple heatmap
#corvals4 <- dba.plotHeatmap(dat.rp2)

#Binding affinity heatmap
#corvals4 <- dba.plotHeatmap(dat.rp2, contrast=1, correlations=FALSE)

Session information

R version 3.5.0 (2018-04-23)

**Platform:** x86_64-apple-darwin15.6.0 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: parallel, stats4, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: bindrcpp(v.0.2.2), DT(v.0.4), dplyr(v.0.7.6), plyr(v.1.8.4), pander(v.0.6.1), DiffBind(v.2.8.0), SummarizedExperiment(v.1.10.1), DelayedArray(v.0.6.1), BiocParallel(v.1.14.1), matrixStats(v.0.53.1), Biobase(v.2.40.0), GenomicRanges(v.1.32.3), GenomeInfoDb(v.1.16.0), IRanges(v.2.14.10), S4Vectors(v.0.18.3) and BiocGenerics(v.0.26.0)

loaded via a namespace (and not attached): amap(v.0.8-16), colorspace(v.1.3-2), rjson(v.0.2.20), hwriter(v.1.3.2), rprojroot(v.1.3-2), htmlTable(v.1.12), XVector(v.0.20.0), base64enc(v.0.1-3), rstudioapi(v.0.7), ggrepel(v.0.8.0), bit64(v.0.9-7), AnnotationDbi(v.1.42.1), splines(v.3.5.0), geneplotter(v.1.58.0), knitr(v.1.20), Formula(v.1.2-3), Rsamtools(v.1.32.0), annotate(v.1.58.0), cluster(v.2.0.7-1), GO.db(v.3.6.0), pheatmap(v.1.0.10), graph(v.1.58.0), compiler(v.3.5.0), httr(v.1.3.1), GOstats(v.2.46.0), backports(v.1.1.2), assertthat(v.0.2.0), Matrix(v.1.2-14), lazyeval(v.0.2.1), limma(v.3.36.2), acepack(v.1.4.1), htmltools(v.0.3.6), prettyunits(v.1.0.2), tools(v.3.5.0), gtable(v.0.2.0), glue(v.1.2.0), GenomeInfoDbData(v.1.1.0), Category(v.2.46.0), systemPipeR(v.1.14.0), ShortRead(v.1.38.0), Rcpp(v.0.12.17), Biostrings(v.2.48.0), gdata(v.2.18.0), rtracklayer(v.1.40.3), stringr(v.1.3.1), gtools(v.3.8.1), XML(v.3.98-1.11), edgeR(v.3.22.3), zlibbioc(v.1.26.0), scales(v.0.5.0), hms(v.0.4.2), RBGL(v.1.56.0), RColorBrewer(v.1.1-2), BBmisc(v.1.11), yaml(v.2.1.19), gridExtra(v.2.3), memoise(v.1.1.0), ggplot2(v.2.2.1), rpart(v.4.1-13), biomaRt(v.2.36.1), latticeExtra(v.0.6-28), stringi(v.1.2.3), RSQLite(v.2.1.1), genefilter(v.1.62.0), checkmate(v.1.8.5), GenomicFeatures(v.1.32.0), caTools(v.1.17.1), rlang(v.0.2.1), pkgconfig(v.2.0.1), BatchJobs(v.1.7), bitops(v.1.0-6), evaluate(v.0.10.1), lattice(v.0.20-35), purrr(v.0.2.5), bindr(v.0.1.1), labeling(v.0.3), GenomicAlignments(v.1.16.0), htmlwidgets(v.1.2), bit(v.1.1-14), tidyselect(v.0.2.4), GSEABase(v.1.42.0), AnnotationForge(v.1.22.0), magrittr(v.1.5), sendmailR(v.1.2-1), DESeq2(v.1.20.0), R6(v.2.2.2), gplots(v.3.0.1), Hmisc(v.4.1-1), DBI(v.1.0.0), foreign(v.0.8-70), pillar(v.1.2.3), nnet(v.7.3-12), survival(v.2.42-4), RCurl(v.1.95-4.10), tibble(v.1.4.2), crayon(v.1.3.4), KernSmooth(v.2.23-15), rmarkdown(v.1.10), progress(v.1.2.0), locfit(v.1.5-9.1), grid(v.3.5.0), data.table(v.1.11.4), blob(v.1.1.1), Rgraphviz(v.2.24.0), digest(v.0.6.15), xtable(v.1.8-2), brew(v.1.0-6) and munsell(v.0.5.0)